    % Matlab code to estimate a flexible mixed logit model in wtp space.
% Written by Kenneth Train, first version Oct 21, 2015, 
% Heather Snell added discounting component April 29, 2019.
 format shortG
 clear all

% Do not change the next line. It sets the global variables.
global NP XMAT PROBS Z NZ NCS NDRAWS NV BETAS ZTYPE CrossCorr YesGPU COEF NROWS discountfun T modelname subset
% OUTPUT FILE
discountfun = "QuasiHyper";
subset = 0; %for full data set
modelname = 'QH';
T=65; %max time periods;
diary off
delete QH.out
diary QH.out
NKnots=8;
filename = 'QH.workspace.mat';
CrossCorr=1;
NBins=20;
WantHessian=1;

WantBoot=0;
NReps=250;
NDRAWS=1000;
NGridPts=6001;

% ZTYPE=1 for polynomials.
%      =2 for step functions
%      =3 for splines
ZTYPE=3;

PolyOrder=8;
NLevels=6;  


if discountfun=="Harvey" % 
R_Range=[0.75 2.3542];
P_Range=[0 4.25];
WTP_Range=[-8.5584    7.1060;
   -7.1937    9.1447;
   -0.6696    1.0361;
   -2.1644    1.9384;
    0.1814    0.2129;
   -0.3655    0.9043;
   -0.0723    0.2217];

elseif discountfun=="HM" % 
    R_Range=[1.32 2.1];
    P_Range = [0 3.36];
    WTP_Range=[-9  10;%cap
               -8.7103  10.602;%sw
               -0.43051  0.7104;%buffer
               -1.5578  1.4062;%infra
                0.11959  0.20079;%jobs
               -0.33872  1.0456;%quality
               -1.359  1.5378];%wlife  
elseif discountfun=="Exponential" % 
    R_Range = [0.083333, 1.1667];
    P_Range = [0, 2];
    WTP_Range=[-10 12;
                -7.2819 10.576;
                -2.6147 4.5;
                -13.562 9.7886;
                  0.000 0.672;
                -11.958 13.952;
                -12.797 16.486];

elseif discountfun=="QuasiHyper" %
    R_Range = [0.15, 1];
    P_Range = [0.2, 1]; %actually the beta range 
    WTP_Range=[ 0  0.25; %price range
               -12	 10;
                -16	 16;
                -3	 6;
                -10.0187 10;
                 0.6	 0.7174;
                -5.7674	 7.75;
                -3.0512	 4.2668];
               

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Below should not change %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(subset==0)
    load('discounting.data.mat');  %This loads XMAT with the variables below
elseif subset==1
    XMAT = csvread('subset1.csv',1,0);
elseif subset==2
    XMAT = csvread('subset2.csv',1,0);
elseif subset==3
    XMAT = csvread('subset3.csv',1,0);
elseif subset==4
    XMAT = csvread('subset4.csv',1,0);
end


FirstRun = 1; %This is used and changed in Bootstapping, but don't change here
Nalts = 3; %number of alts per question
Nquest = 5; %number of questions each person answers
NROWS = size(XMAT, 1);
NP = NROWS/(Nalts*Nquest);
NCS = NP * Nquest; 

% 1. Person number (1-NP)            MUST BE THIS. DO NOT CHANGE.
% 2. Choice situation number (1-NCS) MUST BE THIS. DO NOT CHANGE.
% 3. Chosen alternative (1/0)        MUST BE THIS. DO NOT CHANGE.
% 4. Price 
% 5. CAP
% 6. SW
% 7. Buffer
% 8. Quality
% 9. Jobs
% 10. Infra
% 11. Wlife
% 12. treatment
if discountfun=="QuasiHyper"
    NAMES = {'Rate';  'Beta'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'Buffer';
                'Infra'; 'Jobs'; 'Quality'; 'Wlife'};
else
    NAMES = {'Rate'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'Buffer';
    'Infra'; 'Jobs'; 'Quality'; 'Wlife'};
end

IDPRICE=4;

NV = size(NAMES, 1);
%Do not change the next line
COEF = [R_Range;P_Range;WTP_Range];

ThisSeed=1234;
rng(ThisSeed);

%Do not change the following lines. They calculate the number of Z variables
if ZTYPE==1; 
    NZ=PolyOrder*NV;
    if CrossCorr==1;
        if discountfun=="QuasiHyper"
            NZ=NZ+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
        else 
            NZ=NZ+(NV-2)*(NV-3)/2; %will not correlate price or rate
        end
    end
elseif ZTYPE==2
    NZ=(NLevels-1)*NV;
    if CrossCorr==1;
            if discountfun=="QuasiHyper"
               NZ=NZ+2*(NV-3)+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
            else 
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
elseif ZTYPE==3
    NZ=(NKnots+1)*NV;
    if CrossCorr==1;
            if discountfun=="QuasiHyper"
               NZ=NZ+2*(NV-3)+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
            else 
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
end

%Set the starting values for the coefficients of the Z variables
beta=zeros(NZ,1);  

YesGPU=1;

MAXITERS=5000;

%CONVERENCE TOLERANCE.
% Set =[] to use matlab defaults. default is 1e-6
PARAMTOL=[];
LLTOL=[];
tic;
if discountfun ~= "QuasiHyper"
    CreateDrawsDisWtpProbs;
    CheckRate;
    if check_ok==1
       CreateZ;
       disp(' ');
       disp(['Data setup took ' num2str(toc./60) ' minutes.']);
       disp(' ');
       DoEstimationRate;
        paramhat_orig = paramhat;
        Figure_Output9;
        if WantBoot==1
            Boot
        end
    end
else
    CreateDrawsDisWtpProbsQH;
    CheckRate;
    if check_ok==1
        CreateZQH;
        disp(' ');
        disp(['Data setup took ' num2str(toc./60) ' minutes.']);
        disp(' ');
        DoEstimationRate;
        paramhat_orig = paramhat;
        Figure_Output10;
        if WantBoot==1
           Boot
        end
    end

end

save(filename)
diary off
